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Abstract 

An efficient method to account for the chemi- 
cally frozen thermodynamic and transport properties 
of air in three dimensional Navier-Stokes calculations 
has been demonstrated. This approach uses an explic- 
itly specified equation of state (EOS) so that the fluid 
pressure, temperature and transport properties are di- 
rectly related to the flow variables. Since the pressure is 
explicitly known as a general function of the flow vari- 
ables no assumptions are made regarding the pressure 
derivatives in the construction of the flux Jacobians. 
The method is efficient since no sub-iterations are re- 
quired to deduce the pressure and temperature from 
the flux variables and allows different equations of state 
to be easily supplied to the code. The flexibility of the 
EOS approach is demonstrated by implementing a high 
order TVD upwinding scheme based upon flux differ- 
encing and Van Leer’s flux vector splitting. The EOS 
approach is demonstrated by computing the hypersonic 
flow through the corner region of two mutually perpen- 
dicular flat plates and through a simplified model of a 
scramjet module gap-seal configuration. 

I. Introduction 

The capability to compute chemically frozen sin- 
gle stream flows using an explicitly specified equation of 
state (hereafter Referred to as EOS) has been added to 
the PARC3D Navier-Stokes solver. The PARC codes 
originated as the AIR2D/3D solvers by Pulliam and 
Steger (Ref. 1), which used the ADI algorithm pre- 
sented by Beam and Warming (Ref. 2). To increase 
the computational efficiency and to reduce storage re- 
quirements, Pulliam (Ref. 3) diagonalized the block 
implicit operators of the original Beam- Warming ADI 
algorithm, resulting in the ARC2D/3D codes. Cooper 
(Ref. 4) further modified the ARC codes to be used 
in a production environment for solving the propulsion 
oriented applications at the Arnold Engineering Devel- 
opment Center, resulting in the PARC codes. Coirier 
(Ref. 5) replaced the implicit operator with the 
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LU-SGS (Lower-Upper Symmetric Gauss-Seidel) algo- 
rithm and demonstrated this scheme for hypersonic 
corner and gap-seal problems. This scheme is based 
upon an LU decomposition method proposed by Jame- 
son and Turkel (Ref. 6) and demonstrated by Jameson 
and Yoon for the Euler equations (Ref. 7). In addi- 
tion, Shuen and Yoon (Ref. 8) and Yu, Tsai and Shuen 
(Ref. 9) have shown the efficiency and robustness of 
this algorithm for solving the Navier- Stokes equations 
fully coupled to the chemical transport equations for 
computing finite rate chemically reacting flow fields. 

For flows where the dissociation and recombination 
rates are much slower than the fluid residence time (i.e.- 
in the limit as the Damkohler number approaches zero) 
the fluid can be modelled as a mixture of inert (chemi- 
cally frozen) gases. When the converse is true (i.e.- the 
Damkohler number approaches infinity) the fluid can 
be modelled as if it were in chemical equilibrium. If 
either of these limitations holds for the flow field to be 
considered, great simplifications can be made and the 
resulting computational work reduced. The present re- 
search adds the capability to the PARC codes to model 
general non-mixing chemically frozen flows using an ef- 
ficient method based upon a generalized equation of 
state. The flexibility of this approach is demonstrated 
by implementing a high order TVD upwinding scheme 
based upon flux differencing (Ref. 10) and Van Leer’s 
flux vector splitting approach (Ref. 11). An equation of 
state modelling chemically frozen air is then generated 
and used in the computations of two different types of 
three dimensional, hypersonic flows: flow through the 
corner region of two mutually perpendicular flat plates 
and flow through a simulated scramjet module gap-seal 
configuration. 

II. Equation of State Formulation 

To account for variable thermodynamic and trans- 
port properties it is necessary to relate the pressure, 
temperature and transport properties to the flow vari- 
ables. In the majority of flow solvers the fluid specific 
heats are assumed to be constant, resulting in a simple 
equation of state, namely: 
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The speed of sound is then easily found as 
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For high enthalpy flows the assumption of constant 
7 can lead to over predictions of the temperatures and 
pressures near stagnation conditions. Real gas effects 
for compressible flow calculations have been made by 
many researchers for mixing and non- mixing flows that 
are considered to be chemically frozen or in chemical 
equilibrium or non-equilibrium (e.g.- Refe. 8,9,12 and 
13). For the frozen and finite rate approaches, once the 
chemical composition is known (by either solving the 
species transport equations with or without the species 
production terms), the pressure and temperature are 
found in much the same manner. The mixture enthalpy 
is formed based upon the chemical composition and the 
individual species enthalpy as 


N 

h(T) = J2 Y *hi ( 6 ) 

t=l 

where hi(T) is typically in the form of polynomials in 
temperature. This is then related to the internal energy 
(which is readily obtained from the flux variables) by 
assuming the fluid is ideal, so that 

e = h(T) - P - (7) 

P 

Then the temperature is found by solving (7) for T 
(typically by Newton’s iteration or an interval halving 
method) and the pressure is found from (2). This ap- 
proach allows mixing or non-mixing chemically frozen 
flows to be calculated in a finite rate chemistry code 
by simply switching off the chemical source terms or 
the species equations completely, but loses efficiency 
by having to perform the temperature inversion sub- 
iterations at each grid point for each global iteration 
level. Not only is more work required to perform these 
sub-iterations, the existence of subroutine calls in the 
do loop structure can devectorize computationally in- 
tensive portions of the code (e.g.- pressure terms in 
the inviscid flux vectors for the residual calculations). 


By having an explicit relationship between the pres- 
sure, temperature and flux variables, significant sav- 
ings in work can be achieved by eliminating these sub- 
iterations and devectorizing subroutine calls. 

Recent work by Shuen, et.al. (Ref. 14) and Liou, 
et. al (Ref. 15) in developing flux splitting methods 
for real gases has made use of an explicitly specified 
equation of state (EOS). This concept was successfully 
applied by Liou (Ref. 16) for real gas calculations of 
a corner flow field using a parabolized Navier- Stokes 
code. By having an explicitly specified EOS, advan- 
tages in flexiblity and efficiency for frozen and equilib- 
rium flows are apparent. No ’’effective 7” approxima- 
tions need to be made in forming flux Jacobians and no 
sub-iteration procedure is needed to find the pressure 
and temperature. In addition, the real gas flux split- 
ting methods in Refs. 15 and 16 can be extended to 
the code once the generalized EOS is operational. 

The EOS’s generated for this work use thermo- 
dynamic and transport data from the NASA Lewis 
Chemical Equilibrium program, CEC (Ref. 17). CEC 
is widely used and is accepted as a valuable tool in 
computing general equilibrium compositions and their 
properties. This program accesses an extensive data 
bank that contains thermodynamic data (in the form of 
polynomials in temperature for individual species spe- 
cific heats, enthalpy and Gibb’s energy) for 741 gaseous 
species, and transport property data (in the form 
of polynomials in temperature for individual species 
molecular viscosities and thermal conductivities) for 
154 gaseous species. For mixtures in chemical equi- 
librium it uses a Gibb’s energy minimization procedure 
that is fully outlined in Ref. 17. 

The chemically frozen EOS’s generated in this 
study are made using a least squares surface fit that 
minimizes the relative error between the fit function 
and the supplied data. Given the chemical composi- 
tion, Yi, the internal energy and density are computed 
over a wide range of pressures and temperatures using 
the thermodynamic data in Reference 17. Then p(/?, e) 
and T(p,e) are found by applying the surface fitting 
procedure. Individual species transport properties are 
calculated from Ref. 18, from which the mixture prop- 
erties are found using the Wilke mixing rule. The fits 
for the mixture laminar viscosity and thermal conduc- 
tivities are then made from this data set. 

Now that p(p, e) and T(p, e) are known, the deriva- 
tives and || can be easily found. These are used in 
the formulation of the flux Jacobians in the LU operator 
and are used in (4) to calculate the speed of sound. By 
explicitly knowing the pressure derivatives no approxi- 
mations need to be made in forming the flux Jacobians, 
such as an ’’effective gamma” approach. In addition, 
when forming the inviscid flux vectors the pressure is 
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immediately available from an evaluation of the equa- 
tion of state. The EOS fits are formulated as state- 
ment functions used by the code and are contained in 
FORTRAN include files. By having the EOS’s as state- 
ment functions, the original vectorized structure of the 
code is unaffected and no sub-iterations are required 
to obtain the pressure and temperature. This allows a 
previously vectorized code to maintain the same basic 
structure and computational efficiency. 

II, a: Chemically Frozen Air Equation of State 


The procedure outlined in the previous section has 
been used to generate an equation of state for air un- 
der the assumption that the fluid is ideal and its com- 
position remains chemically frozen. The mole fraction 
composition of air is taken to be Xjv 2 = 0.78847 and 
Xo 2 ~ 0.21153, which is used in conjunction with the 
thermodynamic data from CEC to generate the den- 
sity and internal energy over an extremely wide pres- 
sure and temperature range. The CEC data for en- 
thalpy used to generate the pressure/temperature map 
has been calculated by statistical mechanics over the 
temperature range of 200 < T(Kelvin) < 15000 (Ref. 
19). Since the surface fitting procedure used in the 
generation of the EOS will have the largest errors near 
the boundaries of the domain, an extremely wide pres- 
sure and temperature range was used to compute the 
density and internal energy. The density and internal 
energy were computed given the pressure and temper- 
ature over the range of 1.0A10 -4 < P(Atm) < 100.0 
and 20 < T(Kelvin) < 15000. (For T < 200(Kelvin) 
the specific heats are assumed to be equal to their val- 
ues at T = 200. Great care was taken in choosing the 
surface fitting base functions to keep a high level of fit- 
ting accuracy over the entire fitting range. The root 
mean square of the relative error for the pressure, tem- 
perature, laminar viscosity and thermal conductivity 
fits was (respectively) 0.42, 0.34, 1.7 and 1.2 percent. 
Appendix B contains the surface fitted base functions 
and coefficients for the fits of pressure, temperature, 
laminar viscosity and thermal conductivities. For com- 
paritive purposes, an EOS for calorically perfect air was 
also generated following (1), (2) and (5) using y = 1.4. 

Figures 1 and 2 compare the pressures and tem- 
peratures for the calorically perfect and thermally per- 
fect EOS’s used in the study over a broad temperature 
and pressure range. Figure 3 compares the laminar vis- 
cosities obtained from Sutherland’s law, data from Ref. 
20 and the present data. Figure 4 compares the ther- 
mal conductivities assuming constant Prandtl number 
(P r = 0.72), C p and Sutherland’s law for the viscos- 
ity, data from Ref. 20 and the present data. Figure 1 


shows the that the differences in pressure between the 
calorically perfect air and chemically frozen air EOS’s 
are small, but not minor, while Figure 2 shows that the 
temperature deviations can be quite significant in the 
higher internal energy ranges. Figures 3 and 4 show 
that the fitted data follows the NBS standards data 
more closely than Sutherland’s law for both the viscos- 
ity and the thermal conductivity. 

Ill, Numerical Solution Procedure 

The flow solver used in this study solves the 
Reynolds averaged Navier-Stokes equations cast in di- 
vergence form evaluated in a generalized curvilinear co- 
ordinate system. The equations are solved implicitly 
using the LU-SGS algorithm for approximate Newton 
iteration and are discretized in a finite difference for- 
mulation. The original discretization of the convective 
terms is central differencing with added 1st and 3rd dif- 
ference dissipation in the manner due to Jameson (Ref. 
21). The coefficients controlling the levels of dissipation 
are kept at their default levels (/c 2 = 0.25 and /c 4 = .01) 
for all of the computations performed in this research. 
To aid in the shock capturing ability of the code, a sec- 
ond order TVD upwinding scheme based upon flux dif- 
ferencing (Ref. 10) and Van Leer’s flux vector splitting 
(Ref. 11) has been formulated using the generalized 
equation of state formulation. The following sections 
briefly describe the implicit and explicit portions of the 
solution scheme and show how the generalized equation 
of state is implemented. 

IILa: LU-SGS Algorithm 

The LU-SGS (Lower-Upper Symmetric Gauss - 
Seidel) scheme was first demonstrated by Jameson and 
Yoon (References 22 and 23) in solving the Euler and 
Navier Stokes equations. For computing flows with fi- 
nite rate chemistry this scheme has been shown to be 
very efficient by implicitly solving the Navier-Stokes 
equations fully coupled to the chemical transport equa- 
tions (References 8,9 and 24). For solving the Euler 
or Navier-Stokes equations alone, the resulting scheme 
is scalar diagonal, and requires only two sweeps with 
scalar inversions for both two and three dimensional 
problems. When source terms are included implicitly, 
the approach results in block diagonal operators. Fol- 
lowing the procedure outlined in Reference 5 the hybrid 
LU-SGS scheme can be derived from the following. 

A prototype implicit scheme for solving the Navier- 
Stokes equations in delta form may be written as 

[/ + At(D{A A D V B + D c C)]Aq = - AtR (8) 


3 


where and are difference operators that ap- 

proximate d^drj and d ^ and A, B and C are the inviscid 
flux Jacobians 


a _ 9E dF d_G 
A ~ dq ,B ~ dq ,C ~ dq 


( 9 ) 


If the flux Jacobians in (8) are approximated so that 
they can be split into and matrices, where the 
eigenvalues of the ”+” matrices are non-negative and 
the eigenvalues of the matrices are non-positive, 
and then differenced according to their sign, (8) may 
be written as 


[I + At(V ( A+ +A(A~ +V n B + + A n B~ (10) 

+V ( C + + A ( C")] = —AtR 

The choice of conditioning used to form the ap- 
proximated Jacobians is very important, and follow- 
ing the construction used in Reference 25, the approx- 
imated Jacobians are formed as 


A ± _ A±r A I ± = B±r B I ± = C±T c[ 
2 2 2 

where 


rA,B,c = 0™ax(\\ A>B ,c\) (12) 

where /? > 1.0 and A a,b,c are the eigenvalues of the flux 
jacobians. Therefore, using the conditioning in (11), 
the unfactored scheme may be written as 


[7(l+A<(r yl + rB + r c ))+A<( J 4- +1 + Bt +1 +C | - 1 (13) 

-Af-i - B +_ , - C+JjAg = —AtR 

where (j,k,l) refers to the grid point index. (When an 
index is omitted it is assumed to be at j,k or 1). The 
hybrid LU- SGS scheme is then found by factoring (13) 
into L,D and U operators so that 


LD~ 1 UAq = -AtR 

(14) 

L = I[l + A t{r A + r B + rc)] 

(15) 

-A <(A+_ 1 + B+_ 1 + C+ 1 ) 


D~ l = 7[1 + At(r A + r B + r c )] 

(16) 

U = 7[1+ A t(r A + r B + rc)] 

(17) 

+A t(Aj +1 + £^ +1 + Cf + J 



All of the calculations performed here are only con- 
cerned with obtaining the steady state, so an approxi- 
mate Newton iteration formulation is formed by divid- 
ing (13) by At and taking the limit as At approaches 
oo, resulting in 

LUAq = -(r A + r B + r c )R (18) 

L = [7(r„ + r B + r c ) - Af_ , - B+_ , - C+J (19) 
U = [/(r^ + r B + r c ) + AJ +1 + B^ +1 + (20) 

The solution of (18) is fully vectorizable on planes 
where j+k+l= const ant and is efficient since it requires 
only scalar diagonal inversions. The approach is to first 
solve for an intermediate vector using (19) by sweeping 
in the direction of planes of increasing j-f k+1 and then 
solve for A q using (20) by sweeping in the opposite 
direction. One advantage of the approximate Newton 
iteration formulation is that no choice of time stepping 
or Courant number is required. Results in Reference 5 
show that the scheme is slightly more efficient in pro- 
cessing time per iteration than the diagonalized Beam- 
Warming scheme, which was originally in the code. 

The inclusion of the generalized equation of state 
formulation in the LU-SGS scheme is straightforward. 
Since the pressure derivatives ^ and ^ are directly 
obtainable from the EOS formulation, the inviscid flux 
Jacobians can be found with no ” equivalent gamma” 
assumptions. Therefore, the only modifications to the 
LU-SGS scheme using the generalized equation of state 
formulation is to include the general pressure formula- 
tion and its derivatives in the inviscid flux Jacobians, 
A, B and C in (18), (19) and (20) as well as the speed of 
sound definition in (4). A general inviscid flux vector 
including the pressure derivative terms is included in 
Appendix A. 

IILb: Upwinding Approach 


The upwinding approach used in this research is 
based upon flux differencing (Ref. 10) and Van Leer’s 
flux vector splitting method (Ref. 11). The upwinding 
approach in Reference 10 uses a monotone interpolation 
of the non-normalized nodal point, split inviscid fluxes 
to cell interfaces using a variable order polynomial in- 
terpolation scheme. The numerical fluxes evaluated at 
the cell interfaces consist of the upwinded components, 
and are written as 

* +i = Fij + (21) 

The interpolation used here interpolates cell face nor- 
malized fluxes, simplifying somewhat the approach pre- 
sented in Reference 10. For the ’’positive” flux contri- 


i 



butions, the normalized cell interface flux is monotoni- 
cally interpolated from the nodal values as 


f* +i = Ff + tf 

(22.a) 

= f.;, - 

(22.6) 

where 


6? =C?A~F t + + C}A+F t + 

(23.a) 

67 =Cf A + F7 + C;A~F~ 

(23.6) 

The order of the interpolation scheme and its limiting 
are contained in the coefficients Cf 2 > which are written 
as 


(24 a) 


(24.6) 

cr=#f)i|a -*+§)] 

(24.c) 

+ «-§)] 

(24.d) 

Following Reference 10 the limiters 

in (24) are 

formed based upon the inner products of neighbor- 
ing flux differences. Although this approach does not 
strictly adhere to the traditional method of limiting 
each flux individually, it has been stated to be satisfac- 
tory over a sufficiently wide range of conditions. The 

limiter arguments rf are written as 


+ _ A+FfoA-F? 
r< ~ A- Ffo A- Ff 

(25. a) 

_ _ A+F~ o A~ Fr 
r< ~ A+F~ oA+F~ 

(25.6) 

The limiter function is taken to be the 
modulus (min-mod) function, where 

minimum- 

<f>(r) = maa:(0.0,mm(r, 1.0)) 

(26) 


The variables «, <r and 9 set the order of the 
scheme. For all of the following upwinded calculations, 
2nd order, full upwinding is used so that k = —1.0, 
c 7 ~ 0.0 and 6 = 1.0. Combining the positive and neg- 
ative fluxes, the cell face normalized numerical flux is 
then represented as 

F i+i = + F i+1 ) - - F~) + (<+ - 6~ +l ) 

(27) 


The first term represents a simple numerical average of 
the nodal fluxes while the combination of the first and 
second terms results in a first order upwinding scheme. 
The third terms represent anti-diffusive fluxes which 
increase the order of accuracy of the scheme, and are 
limited to make the scheme monotone. 

After the construction of the numerical fluxes at 
the cell interfaces, the convective contribution to the 
residual from this coordinate direction is then calcu- 
lated as 

Ri = (F6) i+ ,-(F6) i _ i (28) 

where 6 is the cell interface area. The areas and pro- 
jected areas at the cell interfaces are found using the 
transformation metrics evaluated at the nodal points in 
a consistent manner as (for the £ - flux) 



(29. a) 

(U+i = f((;j L ). + (|).+il 

(29.6) 

(«.)<+» = j[(y)< + (y)*.] 

(29.c) 


so that the cell interface area is 

S = yj 5 X 2 + 2 + 6 Z 2 (30) 

Since this scheme operates only on flux differences, 
it is possible to use either a flux vector splitting ap- 
proach or a flux difference splitting approach. Van 
Leer’s flux vector splitting was chosen due to it’s sim- 
plicity and efficiency and for the ease in which the gen- 
eralized equation of state could be included. It should 
be noted that the equation of state formulation does not 
preclude the use of Roe’s flux differencing approach, 
which would include the explicitly specified pressure 
derivatives §£ and in the formulation of the trans- 
formation matrices. The flux vector splitting used here 
is formulated including the pressure and sound speed 
to facilitate using the generalized EOS, so that when 
\M C \ < 1.0 

/ ±^(M C ± l ) 2 \ 

F *(u F2a + «e)) 

F ± = Ft(v- 6 f^(T2a + n e )) (31) 

pH™ f 2° + « c )) 

\ if [tf-m(w c Ta) 2 ] / 

and when 

M c > 1.0 : F + = F, F~ = 0 (32.a) 
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M e < 1.0 : F + = 0, F~ = F (32.6) 

where the contravariaut velocity and Mach number are 
defined as 


inclusion of a high order TVD upwinding scheme are 
motivated by this effort in the hope that a more thor- 
ough understanding of the underlying fluid dynamic 
processes can help the design effort. 


« c = 


6 x u + 6 y v + 6 z w 


a 


(33) 

(34) 


The parameter m in (31) is found by requiring the 
bracketed terms in the energy flux form a perfect 
square, resulting in 


m 


h 


1 + 2 


a 3 u 3 -fv 3 +tu 3 


(35) 


This formulation of m results in an indeterminate ex- 
pression at stagnation points or along no-slip bound- 
aries, so the relation is modified slightly to give 


m — 



(36) 


The inclusion of the generalized equation of state 
formulation into this upwinding scheme is straightfor- 
ward. At each nodal point the flux variables directly 
yield p and e which in turn give p and a 2 directly from 
the equation of state from which the flux vector (31) 
is constructed. This approach using the generalized 
equation of state requires no special sub-iterations or 
procedures to find the pressure, sound speed or static 
enthalpy, and allows the upwinding approach to be fully 
vectorized. The following sections demonstrate the gen- 
eralized equation of state approach for a selected num- 
ber of three dimensional, hypersonic viscous computa- 
tions. 


IV: Corner and Gap Seal Flow Field Study 


Effort is currently underway to investigate the heat 
transfer and aerodynamic loads in and around a scram- 
jet module gap seal configuration. The seal system 
must prevent high temperature and pressure gases from 
leaking through the gaps between the articulating en- 
gine panel walls and stationary splitter panels (Ref. 
26). This system is crucial to the design of hypersonic 
airbreathing propulsion systems and presents a signifi- 
cant challenge to structural designers to build a viable 
seal system. To assist the design effort, the present re- 
search complements the structural analysis underway 
in an attempt to better understand the fluid dynamic 
and heat transfer phenomena in and around the gap 
seal system. The addition of real gas effects and the 


IV.a: Hypersonic Corner Validation 


To obtain some level of confidence in the modelling 
of the gap seal flow the 3D Navier-Stokes solver has 
been validated against experimental data taken from a 
class of flows similar to the gap seal flow (Reference 
5). The gap-seal geometry is similar to two mutu- 
ally perpendicular flat plates (corner) aligned with the 
freestream flow with a ’’groove” situated at the plate 
intersections. Due to this similarity, the 3D code was 
validated against the experimental data of hypersonic 
corner flow taken by Cresci (Reference 27). These tests 
were conducted at a Mach number of 11.8 over a range 
of Reynolds numbers of 0.15 X 10 6 to 0.5 X 10 6 per 
foot. The data presented in Ref. 5 is axial pressures 
and heat transfer rates plotted against the hypersonic 
interaction parameter 


X = 



(37) 


for different lateral locations measured from the corner. 
In addition, surface pressures, heat transfer rates and 
skin friction data comparisons were made as a func- 
tion of lateral distance from the corner at various axial 
locations. This data not only tests the prediction of 
the mutual interaction of the two plate boundary lay- 
ers upon each other but also the self induced leading 
edge pressures caused by the boundary layer growth 
close to the leading edges. 

Comparisons were made against this set of exper- 
imental data using the 2nd order upwinding approach 
outlined here. (The complete set of comparisons using 
the original central difference formulation are available 
in Reference 5. Since the stagnation temperature was 
low ( T stt oo = 944 K), the generalized equation of state 
was not used, so that in essence, this validation effort 
tests the upwinding approach. A relatively coarse grid 
(36 X 41 X 41) was used for plate lengths and widths 
of 8 and 2.5 inches, respectively. The Reynolds num- 
ber per foot was taken to be 0.15 x 10 6 , resulting in 
a Reynolds number based on plate length of 1.0 x 10 5 
from which the flow is assumed to be laminar. The grids 
were clustered normally in the cross- stream planes us- 
ing the Robert’s transformation with 0 = 1.01 and ax- 
ially near the leading edge using 0 — 1.1, with five grid 
points extended ahead of the leading edge. 
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The comparisons to experiment of the upwinded 
computation were nearly identical to the central dif- 
ference results. Since the central difference results al- 
ready computed in Reference 5 were favorable, only a 
selected portion of the comparison to experiment are 
shown here. Figure 5 shows the predicted wall pres- 
sures versus the hypersonic interaction parameter at 
a station 0.125 inches away from the corner intersec- 
tion, while Figure 6 compares the computed and mea- 
sured lateral variation of heat transfer at a location of 
X = 8.3. As can be seen from these figures, the up wind- 
ing approach correctly simulates the mutual interaction 
of the plate boundary layers and the lateral variation 
of the heat transfer to the plates caused by the leading 
edge shock and boundary layer interaction. 

IV. b: Hypersonic Corner Flow 


The hypersonic laminar flow along the intersec- 
tion of two mutually perpendicular flat plates at a 
freestream Mach number of 12 and an altitude of 
100,000 feet was calculated using the calorically perfect 
and thermally perfect air EOS formulations. Both the 
central difference and 2nd order upwinding approaches 
were used. The length of the plates is taken to be 12 
inches while the width of the plates are both taken to 
be 2.5 inches. At this altitude the freestream tempera- 
ture and pressure are 418.79 Rankine and 23.085 psf., 
resulting in a Reynolds number based on plate length of 
1.226 x 10 6 , permitting a laminar calculation. To high- 
light the real gas effects, adiabatic wall boundary condi- 
tions are taken along the plates, while grid line extrap- 
olation is taken at the outflow and lateral boundaries. 
Although the application of adiabatic wall boundary 
conditions results in unrealistically high wall tempera- 
tures, it is an extreme test of the code robustness and 
the installation of the real gas effects. Since the ther- 
mally perfect EOS formulation allows for a variation in 
the fluid specific heats while the calorically perfect for- 
mulation does not, the difference in wall temperatures 
between these formulations should be readily appar- 
ent. The grid is clustered normally to the walls using a 
Robert’s transformation with 0 =1.01 and is clustered 
axially near the leading edge using 0 =1.1. 31 grid 
points are used in the axial direction on the plates with 
5 points extended ahead of the plate leading edges using 
a grid spacing equal to that of the first cell at the plate 
leading edge while 41 points are used in the direction 
normal to each plate. 

Figure 7 illustrates the mutual interactions of 
the plate boundary layers and embedded leading edge 
shock structure and the downstream growth of the in- 
teraction region by showing axial velocity contours in 


planes perpendicular to the plates at different down- 
stream axial locations. Figure 8 shows the predicted 
lateral variation of wall temperature approximately 6 
inches downstream of the leading edges for the calori- 
cally perfect gas equation of state for both the central 
difference and upwinded discretizations, while Figure 
9 shows the chemically frozen air equation of state re- 
sults. The peak wall temperature locations are prac- 
tically identical for the two EOS formulations, with 
the central difference scheme consistently predicting 
lower peaks than the upwinded scheme. These fig- 
ures also show that the wall temperature in the non- 
interactive region computed by the central difference 
scheme is lower than that computed by the upwind 
scheme. The non-interactive region of the corner flow- 
field is far enough downstream of the leading edge so 
that the leading edge effects of the primary plate are 
absent, and far enough laterally from the corner so that 
the opposing plate leading edge shock effects are also 
absent. This difference in computed non-interactive 
wall temperatures indicates that the upwinding scheme 
appears to be less diffusive than the central difference 
scheme. Figure 9 shows the marked reduction in wall 
temperature by using the chemically frozen air equa- 
tion of state, which in essence allows for variation of the 
specific heats. Similar trends between the differencing 
schemes as in the calorically perfect computations are 
also present here. 

Differences between results of the two differencing 
schemes are apparent from a detailed examination of 
cross stream contours of density and cross stream ve- 
locity vector plots. Figure 10 compares contour plots of 
density from both the central difference and upwinded 
discretization schemes using the chemically frozen air 
equation of state. These density contours are taken 
in a plane perpendicular to the free stream at a dis- 
tance approximately 3.5 inches from the leading edge 
of the plates. The central differencing results, shown 
in Figure 10. a. indicate a merging of the leading edge 
shock and the plate boundary layer and an oscillatory 
shock/shock interaction region. The upwinded results, 
in Figure 10. b, indicate a more clearly defined and sep- 
arate leading edge shock and boundary layer, and a 
cleaner shock/shock interaction. These subtle differ- 
ences between the two discretization schemes produce 
rather different results in the cross stream velocity field, 
which is shown in Figure 11. The central difference re- 
sults, in Figure 11. a, show the beginning of a region 
of cross stream separation, which results in the down- 
stream development of a weak vortical flow just inboard 
of the embedded shock. The upwinded results, in Fig- 
ure ll.b, indicate a much weaker cross stream separa- 
tion, resulting in a very weak vortical flow downstream. 
If the upwinding scheme is indeed less dissipative than 
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the central difference scheme, as is indicated in the adi- can see that the basic geometry is similar to two mutu- 

abatic wall temperatures, then it is possible that this ally perpendicular flat plates with a ’’groove” embedded 

vortical flow is induced by the numerical dissipation. in the corner. Studies show that this groove, or recess, 

This finding of numerically induced vortical flow has fills up with low momentum and low energy fluid, which 

been previously investigated in detail for conical solu- reduces the heat loads to the seal. Since this configu- 

tions to the Euler equations for flow over a delta wing ration is intended to be used at high Mach numbers, 

(Ref. 28). The results of that study indicated that the the effects of real gas properties and the inclusion of 

existence and strength of the computed cross stream a high order upwinding scheme may be important in 

vortical flow on the lee side of a blunt tip delta wing characterizing the heat and aero loads upon the seal, 

was sensitive to the level of the 2nd and 4th order dis- Simulations of the gap seal flowfield were made at 

sipation coefficients. This study indicated that with a freestream Mach number of 10 and an altitude of 

lower numerical dissipation, the vortical flow disap- 100,000 feet. The gap seal recess measured 0.25 by 0.25 

pears. Other computational results have indicated the inches in the corner region of the two mutually perpen- 

existence of a vortical structure in corner flow fields, but dicular flat plates. The plates lengths and widths were 

the findings here indicate that the strength is influenced taken to be 10 and 1 inches, respectively. At this alti- 

by the level and type of numerical dissipation. Clearly, tude and freestream Mach number, the Reynolds num- 

more detailed experimental data and more refined nu- ber based on plate length is 8.5139 x 10 5 , permitting 

merical experiments are necessary to confidently pre- a laminar calculation. Grid points were clustered nor- 

dict the evolution and behavior of the corner induced mally to the recess walls and the corner walls and were 

vortical flow. clustered axially towards the leading edge with 5 points 

The relative efficiency using the generalized equa- extended ahead of the leading edge. The calculations 

tion of state can be shown by comparing the process- were performed ona36X51X51 grid, of which ap- 
ing rate per iteration of the different methods. All proximately 78.5 percent of the grid points were active 

of the computations were made on a CRAY-YMP us- in the calculation. As in the previous case, adiabatic 

ing the cft77 compiler. The central difference com- wall boundary conditions were used to highlight the 

putations for the calorically perfect formulation and real gas effects. Similar to the corner case, central dif- 

the thermally perfect formulation took 2.4575 x 10 -5 ferencing and 2nd order upwinding were used for both 

and 2.7672 x 10“ 5 seconds/iteration/grid point, respec- the calorically and thermally perfect equation of state 

tively. This represents a 12.6 percent increase in com- formulations. 

puting time to include the generalized equation of state. Figure 12 shows contours of axial velocity in planes 

The 2nd order upwinding computations for the calor- perpendicular to the freestream at different axial loca- 

ically perfect and thermally perfect formulations took tions of the gap seal configuration. This figure shows 

3.0829 x 10 5 and 3.5453 x 10 5 seconds/iteration/grid the ra pid thickening of the inter-recess boundary layers 

point, respectively. This represents a 15 percent in- due to the over pressures caused by the recess leading 

crease in computing time to include the generalized edge shocks and their subsequent reflections off of the 

equation of state. These processing rates also indicate recess walls. Close inspection of the entrance region 

that the addition of the 2nd order upwinding scheme j n the recess reveals a corner-like flow structure. This 

was made at a cost of approximately 25 percent over similarity to the corner flow only lasts until the inter- 

that of the central differencing scheme. recess leading edge shocks impinge upon the boundary 

layers, causing a rapid deceleration of the flow in the re- 
IV. c: Hypersonic Gap Seal Simulation cess. The mutual interactions of the recess leading edge 

shocks and the expulsion of fluid from the recess area 
results in complicated shock-shock and shock-boundary 
To better understand the heat transfer loads upon layer interactions, 

the gap seal region in a scramjet combustor it is nec- Figures 13 and 14 compare cross-stream pressure 

essary to have an understanding of the flow field in contours computed with the thermally perfect equa- 

and above the seal region. This seal prevents hot, high tion of state for the central difference and upwinding 

pressure gases from leaking through the juncture of the schemes. Figures 13. a and 13.b show pressure con- 

articulating engine panels and the stationary sideplate tours computed by the central difference and upwinding 

walls in a scramjet module (see Ref. 26 for a detailed schemes in a plane located approximately 0.6 inches 

description of the gap seal methodology). Preliminary from the leading edge of the gap seal configuration, 

work has begun to characterize this flowfield and the This figure shows the sideplate leading edge shock (on 

resulting heat loads upon the seal region using Navier- the left), the three recess leading edge shocks interac- 

Stokes technology (Ref. 5). Referring to Figure 12, one tion upon each other and the merging of the cowl in- 


8 


duced leading edge shock and the shock from the right 
recess wall. The most noticeable difference between the 
two schemes is the smearing of the cowl leading edge 
shock to the cowl wall and the oscillations of the pres- 
sure near the convex corner of the recess/cowl junc- 
tion by the central difference approach. Much more 
noticeable differences are apparent downstream from 
this plane, as shown in Figures 14. a and 14.b. At this 
location, the sideplate leading edge shock has traversed 
out of the recess area and merged with the cowl induced 
leading edge shock. This results in essentially two sepa- 
rate shock-shock interaction regions: the cowl and side- 
plate leading edge shock’s interactions (the upper in- 
teraction region), and the recess leading edge shock’s 
interactions (the interaction still in the recess). Figure 
14.a (from the central difference scheme) shows smear- 
ing of the cowl leading edge shock and large oscilla- 
tions that persist near the cowl as well as smaller os- 
cillations inboard of the sideplate induced leading edge 
shock. Figure 14. b (from the upwinding scheme) has 
less smearing of the cowl leading edge shock, predicts 
a more compact cowl/sideplate shock-shock interaction 
and does not have the severe oscillations near the cowl 
wall as the central difference scheme’s results. 

Previous computations using the central difference 
scheme of the gap seal flow field in Reference 5 indicate 
a vortical flow on the cowl plate, situated just outboard 
of the recess. Both the central difference and upwinded 
computations performed here predict a similar vorti- 
cal flow, although the strength of this flow is reduced 
considerably with the upwinded computations. Fig- 
ures 15. a and 15.b illustrate this by showing cross flow 
velocity vectors at a location approximately 6 inches 
from the leading edges of the gap seal configuration. 
The central difference results, in Figure 15. a, show a 
much stronger and larger vortical structure than the 
upwinded results, shown in Figure 15. b. This phenom- 
ena is similar to the results in the corner flow simu- 
lation, in that although both schemes predict a vor- 
tical flow, the central differencing scheme predicted a 
stronger vortical flow. 

Figure 16 compares wall temperatures along the 
centerline of the recess area for the calorically and ther- 
mally perfect equations of state for the upwinding and 
central difference schemes. Since the thermally perfect 
equation of state essentially allows for a variation in the 
fluid specific heats while the calorically perfect equa- 
tion of state does not, the wall temperature along the 
centerline is reduced considerably. As in the previous 
corner case, the relative efficiency using the generalized 
equation of state is shown by comparing the process- 
ing rate per iteration of the different methods. As be- 
fore, the computations were made on a CRAY-YMP 
using the cft77 compiler. The central difference com- 


putations for the calorically perfect formulation and 
the thermally perfect formulation took 2.3078 x 10”' 5 
and 2.5832 x 10” 5 seconds/iteration/grid point, respec- 
tively. This represents a 11.9 percent increase in com- 
puting time to include the generalized equation of state. 
The 2nd order upwinding computations for the calor- 
ically perfect and thermally perfect formulations took 
2.9329 x 10“ 5 and 3.4121 x 10~ 5 seconds/iteration/grid 
point, respectively. This represents a 16.3 percent in- 
crease in computing time to include the generalized 
equation of state. 

IV: Summary and Conclusions 

An efficient method to account for limited real gas 
effects in general Navier-Stokes computations has been 
presented. This method specifies the fluid pressure, 
temperature, laminar viscosity and thermal conductiv- 
ity as explicit functions of the density and internal en- 
ergy, which are simply related to the flux variables. By 
having an explicitly specified equation of state (EOS), 
no sub- iterations are required to deduce the pressure 
and temperature from the flux variables. The EOS ap- 
proach was demonstrated by generating an equation 
of state for air under the assumption that it is ther- 
mally perfect, to allow for the variation of fluid specific 
heats with temperature. The EOS’s were made avail- 
able to the solver as statement functions, which allows 
a previously vectorized code to maintain its same basic 
structure and efficiency. 

The flexibility of the EOS approach was demon- 
strated by including real gas effects in both the implicit 
and explicit portions of the numerical scheme. Since 
derivatives of the pressure with respect to the flux vari- 
ables are immediately available from the EOS, they are 
included in the implicit portion of the scheme without 
any assumptions. The central difference scheme with 
1st and 3rd difference dissipation originally present in 
the code was easily modified with the EOS formula- 
tion. A second order, TVD upwinding scheme was then 
added to the solver, which was also readily modified 
using the EOS approach. The extension to upwinding 
was then favorably validated against experimental data 
from a hypersonic corner flow experiment. 

To test the installation of the real gas effects, and 
to compare the discretization schemes, two different 
flow fields were computed with and without the EOS. 
The first simulation was of the hypersonic flow through 
the corner region of two mutually perpendicular flat 
plates aligned with the freestream. Marked differences 
in wall temperatures were noted between the calori- 
cally perfect and thermally perfect formulations, as well 
as significant differences between the central difference 
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and upwinding schemes. The second simulation was of 
the hypersonic flow in and around a scramjet gap seal 
configuration, which was represented by a corner ge- 
ometry with a recessed seal, or groove, taken out along 
the plate intersections. Major similarities and differ- 
ences between the corner and gap seal flow fields were 
noted, as well as differences in the shock-shock inter- 
action phenomena predicted with the central difference 
and upwinding schemes. Relative efficiencies between 
the calorically perfect and thermally perfect gas for- 
mulations were made by comparing processing rates. 
These comparisons indicated that the thermally per- 
fect EOS was approximately 12 percent slower than the 
calorically perfect gas mode, and that the upwinding 
scheme took approximately 25 percent more time than 
the central differencing scheme. 


Appendix A: Generalized Flux Jacobians 


The Jacobians of the generalized inviscid flux vec- 
tors explicitly including the pressure derivative terms 
from the equation of state are listed below. 

Let the flux variables be 




( p \ 

pu 

pv 

pw 

\E t / 


and the flux vector be 


(A. 1) 
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pu e \ 
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\ (£t + P) u c / 
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(AS) 
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u c = 6 x u -f 6 y v -f 6 z w (A.3) 


Appendix B: Frozen Air Equation of State 


where the metric terms 6 X , 6 y and S 2 take the form for 
the £ flux 



(A A) 


The following shows the equation of state fits gen- 
erated for chemically frozen air using the thermody- 
namic and transport data accessed by the CEC pro- 
gram. The supplied curve fit functions and arguments 
have been non-dimensionalized by the following: 


Then, the flux Jacobian can be formed as 


A = 


dF 

dq 


(AS) 


Poo =101325.0( — r) (5.1. a) 

m z 

^ =m.0(Kelvin) (5.1.6) 
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Poo =1.17196272(^1) (B.l.c) 

m 2 

eoo =8.64571924 x 10 4 (^-) ( B.l.d ) 

Poo =180.16133(^5) (S.l.e) 

A TO =12.619813(£|) (5.1./) 


The EOS functions for pressure, temperature, laminar 
viscosity and thermal conductivity are 

p(p, e) =a x e + a 2 pe + a 3 + a 4 pe 2 + a 5 pe 3 + ( B.2.a ) 

a 6 pe 4 + a 7 p 2 e 3 + a s p 2 e 4 + a g pe 5 + 
aioP 2 e 5 + a n p 2 e 6 + a^pc 6 + a\ 3 p 2 e 7 + a\ 4 pe 7 

=a 2 e + a 4 e 2 + a 5 e 3 + a 6 e 4 + 2a 7 pe 3 + (5.2.6) 

2 a 8 pe 4 + a 9 e 5 + 2ai 0 pe 5 + 2anpe 6 + 
ai 2 e 6 + 2a\ 3 pe 7 4- aj 4 e 7 

dp 

— =a x -f a 2 p + 2 a 4 pe + 3 a 5 pe 2 -f 4 a 6 pe 3 + ( B.2.c ) 

Za 7 p 2 e 2 -f- 4agp 2 e 3 4* 5a 9 pe 4 -f 5ai 0 p 2 e 4 + 

^ a n p 2 4* 6ai2pe 5 + 7ai3/? 2 e 6 + 7a\ 4 pe 6 

T(p, e ) =b\e/p 4* b 2 e 4- b 3 /p 4- b 4 e 2 4- b^e 3 4- bee\B.3.a) 
4* b 7 pe 3 4* b$pe 4 4- 6 9 e 5 4- bi 0 pe 5 4- b n pe 6 + 
b\ 2 e 6 4- bispe 7 4- bi 4 e 7 

dT 

~ = ~~ b\e/ p 2 — 6 3 /p 2 4* &7e 3 4* b&e 4 4- 6ioe 5 (B.3.6) 

4- inc 6 4- bi^e 7 


The EOS coefficients a n , 6 n , c n and are tabu 
lated below. 


n 

an 

bn 

1 

-1.64184011 x 10-° 8 

-4.66744619 x 10“ n 

2 

+4.25555113 x 10~ 01 

+4.11592537 x 10" 01 

3 
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4 

-6.01056836 x 10"° 3 

-4.76679827 x 10“° 3 

5 

+1.16745176 x 10~° 4 
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6 
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7 

+1.68428888 x 10~° 7 

+1.35065409 x 10" 08 

8 

-4.49056339 x 10" 09 
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9 

+6.28302532 x 10“° 9 

+3.84665192 x 10~° 9 

10 

+4.16636684 x 10" 11 

+2.65689892 x 10~ 12 

11 

-1.62288924 x 10~ 13 

-9.73342645 x 10" 15 

12 

-1.66714128 x 10 -11 
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13 

+2.26513998 x 10“ 16 

+1.30418711 x lO" 17 

14 

+1.75501828 x 10" 14 

+9.61090540 x 10~ 15 

n 

c rx 

d n 

1 

-5.65406343 x 10" 04 

+1.31392583 x 10" 01 

2 
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+8.37926287 x 10“° 3 

3 
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4 

-1.21977866 x 10' 01 

-2.58670568 x 10~ 01 

5 

+1.59535520 x 10“° 2 
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6 

-6.81970909 x 10“ 04 

-1.33315276 x 10“° 3 


-fo ~b\/ p 4- b 2 4- 26 4 e 4- 36 5 e 2 + 46 6 e 3 4- {B.Z.c) 

3b 7 pe 2 4- 4 b 8 pe 3 4- 5b 9 e 4 4- 5&i 0 /?e 4 4- 
6 b n pe 5 4- 66 12 e 5 + 7b 13 pe 6 + 7b 14 e 6 


p(T) = Cl + c 2 T* + c 3 T + C4 rt+ 
c 5 t 2 + c 6 rt + C7 t 3 

(BA) 

X(T) =d! + d 2 T$ + d 3 T + d 4 rf + 
d 5 T 2 + d 6 T* + d 7 T 3 

(5.5) 
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Figure 1. Equation of State Comparison: 
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Figure 2. Equation of State Comparison: 
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Figure 3. Equation of State Comparison: 
Laminar Viscosity 


Figure 4. Equation of State Comparison: 
Thermal Conductivity 
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Figure 15.a Cross Flow Velocity Vectors, 

Central Differencing: Gap Seal Case 


Figure 15.b Cross Flow Velocity Vectors, 

Flux Differencing: Gap Seal Case 
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Figure 16. Recess Centerline Wall Temperatures 
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sub-iterations are required to deduce the pressure and temperature from the flux variables and allows different 
equations of state to be easily supplied to the code. The flexibility of the EOS approach is demonstrated by 
implementing a high order TVD upwinding scheme based upon flux differencing and Van Leer’s flux vector 
splitting. The EOS approach is demonstrated by computing the hypersonic flow through the corner region of two 
mutually perpendicular flat plates and through a simplified model of a scramjet module gap-seal configuration. 
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